* Independently generated unit root outcome and white noise predictors DGP

// Basics
cscript
cls
set matsize 11000
parallel setclusters 8

capture program drop NG
program define NG, rclass

// Monte Carlo Basics
version 12.1
drop _all
syntax [, obs(integer 1) ]
set obs `obs'


// Time Series Basics
egen time = fill(0 1 2)
tsset time

// Data Generating Process

* Stationary predictor
gen x1 = rnormal()

* Stationary predictor
gen x2 = rnormal()

* Generate outcome
gen y = rnormal()
replace y = L.y + rnormal() if time > 0

// Model Estimation

local phi_y = 1
* Estimate GECM //correct
reg D.y L.y D.x1 L.x1 D.x2 L.x2
local gecm_true_phi = `phi_y' - 1
scalar phi_gecm_test = (_b[L1.y] - `gecm_true_phi') / _se[L1.y]
return scalar phi_gecm_t1 = 2*ttail(e(df_r),abs(phi_gecm_test)) < 0.05
test D.x1 D.x2
return scalar type1_x_gecm = r(p) < 0.05
test L.x1 L.x2
return scalar type1_lx_gecm = r(p) < 0.05
test D.x1 L.x1
return scalar type1_x1_gecm = r(p) < 0.05
test D.x2 L.x2
return scalar type1_x2_gecm = r(p) < 0.05
scalar phi_gecm_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_gecm_power = 2*ttail(e(df_r),abs(phi_gecm_test_power)) < 0.05

* Estimate ADL //correct
reg y L.y x1 L.x1 x2 L.x2
scalar phi_adl_test = (_b[L1.y] - `phi_y') / _se[L1.y]
return scalar phi_adl_t1 = 2*ttail(e(df_r),abs(phi_adl_test)) < 0.05
test x1 x2
return scalar type1_x_adl11 = r(p) < 0.05
test L.x1 L.x2
return scalar type1_lx_adl11 = r(p) < 0.05
test x1 L.x1
return scalar type1_x1_adl11 = r(p) < 0.05
test x2 L.x2
return scalar type1_x2_adl11 = r(p) < 0.05
scalar phi_adl_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi_adl_power = 2*ttail(e(df_r),abs(phi_adl_test_power)) < 0.05


* Estimate ADL (2,2) //correct
reg y L.y L2.y x1 L.x1 L2.x1 x2 L.x2 L2.x2
scalar phi1_adl22_test = (_b[L1.y] - `phi_y') / _se[L1.y]
return scalar phi1_adl22_t1 = 2*ttail(e(df_r),abs(phi1_adl22_test)) < 0.05
test x1 x2
return scalar type1_x_adl22 = r(p) < 0.05
test L.x1 L.x2
return scalar type1_lx_adl22 = r(p) < 0.05
test x1 L.x1
return scalar type1_x1_adl22 = r(p) < 0.05
test x2 L.x2
return scalar type1_x2_adl22 = r(p) < 0.05
scalar phi1_adl22_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi1_adl22_power = 2*ttail(e(df_r),abs(phi1_adl22_test_power)) < 0.05

* Estimate ADL (3,3) //correct
reg y L.y L2.y L3.y x1 L.x1 L2.x1 L3.x1 x2 L.x2 L2.x2 L3.x2
scalar phi1_adl33_test = (_b[L1.y] - `phi_y') / _se[L1.y]
return scalar phi1_adl33_t1 = 2*ttail(e(df_r),abs(phi1_adl33_test)) < 0.05
test x1 x2
return scalar type1_x_adl33 = r(p) < 0.05
test L.x1 L.x2
return scalar type1_lx_adl33 = r(p) < 0.05
test x1 L.x1
return scalar type1_x1_adl33 = r(p) < 0.05
test x2 L.x2
return scalar type1_x2_adl33 = r(p) < 0.05
scalar phi1_adl33_test_power = (_b[L1.y] - 0) / _se[L1.y]
return scalar phi1_adl33_power = 2*ttail(e(df_r),abs(phi1_adl33_test_power)) < 0.05

end
